Mouse analysis

packages <- c("CIMseq", "CIMseq.data", "tidyverse", "scSeqTools")
purrr::walk(packages, library, character.only = TRUE)
rm(packages)

##DATA
load('../data/CIMseqData.rda')
load('../data/sObj.rda')

#rename classes
renameClasses <- function(class) {
  case_when(
    class == "0" ~ "C.Stem.proximal",
    class == "1" ~ "C.Stem.distal",
    class == "2" ~ "SI.Stem",
    class == "3" ~ "C.Goblet.distal",
    class == "4" ~ "SI.Enterocytes",
    class == "5" ~ "SI.Goblet",
    class == "6" ~ "C.Goblet.proximal",
    class == "7" ~ "Enteroendocrine",
    class == "8" ~ "Tufft",
    class == "9" ~ "SI.Paneth",
    class == "10" ~ "Blood",
    TRUE ~ "error"
  )
}

cOrder <- c(
  "C.Stem.proximal", "C.Goblet.proximal", "C.Stem.distal", "C.Goblet.distal",
  "SI.Goblet", "SI.Paneth", "SI.Stem", "SI.Enterocytes",
  "Enteroendocrine", "Tufft", "Blood"
)

getData(cObjSng, "classification") <- renameClasses(getData(cObjSng, "classification"))
features <- getData(cObjMul, "features")
names(features) <- renameClasses(names(features))
getData(cObjMul, "features") <- features
fractions <- getData(sObj, "fractions")
colnames(fractions) <- renameClasses(colnames(fractions))
sObj@fractions <- fractions

plotUnsupervisedClass(cObjSng, cObjMul)

plotUnsupervisedMarkers(
  cObjSng, cObjMul,
  c("Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1"),
  pal = RColorBrewer::brewer.pal(8, "Set1")
)

plotUnsupervisedMarkers(
  cObjSng, cObjMul, "Hoxb13",
  pal = RColorBrewer::brewer.pal(8, "Set1")
)

plotUnsupervisedMarkers(
  cObjSng, cObjMul, "Mki67",
  pal = RColorBrewer::brewer.pal(8, "Set1")
)

classHeatmap(
  data = data.frame(
    gene = rownames(getData(cObjMul, "counts"))[getData(cObjMul, "features")],
    class = names(getData(cObjMul, "features"))
  ),
  counts.log = getData(cObjSng, "counts.log"),
  classes = tibble(
    sample = colnames(getData(cObjSng, "counts")), 
    class = getData(cObjSng, "classification")
  )
)

violinPlot <- function(cpm, genes, classes, class) {
  cpm[genes, ] %>%
    t() %>% 
    matrix_to_tibble("sample") %>%
    gather(gene, cpm, -sample) %>%
    full_join(tibble(sample = colnames(cpm), class = classes), by = "sample") %>%
    ggplot() +
    geom_jitter(aes(class, cpm), height = 0, size = 0.1, alpha = 0.25) +
    facet_wrap(~gene, scales = "free_y") +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90)) +
    labs(title = paste0(class, " genes"))
}


getGenes <- function(
  cpm, classes, features, class, ntop = 10
) {
  fc <- foldChangePerClass(
    cpm[features[names(features) == class], ], 
    tibble(sample = colnames(cpm), class = classes)
  )
  rownames(fc[order(fc[, class], decreasing = TRUE), ])[1:ntop]
}

c <- "C.Stem.distal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

c <- "C.Stem.proximal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

c <- "C.Goblet.distal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

c <- "C.Goblet.proximal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

c <- "SI.Stem"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

c <- "SI.Goblet"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

c <- "SI.Enterocytes"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

c <- "SI.Paneth"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)

Deconvolution

plotSwarmCircos(sObj, cObjSng, cObjMul, weightCut = 0, classOrder = cOrder)

Filtering

Only detected duplicates and triplicates.
Only ERCC estimated cell number <= 4.
Weight cutoff = 5.

adj <- adjustFractions(cObjSng, cObjMul, sObj, binary = TRUE)
samples <- rownames(adj)
rs <- rowSums(adj)
keep <- rs == 2 | rs == 3

plotSwarmCircos(
  filterSwarm(sObj, keep), cObjSng, cObjMul, weightCut = 5, 
  classOrder = cOrder, theoretical.max = 4
)